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We introduce a new quantity to probe the glass transition. This quantity is a linear general- 
ized compressibility which depends solely on the positions of the particles. We have performed a 
molecular dynamics simulation on a glass forming liquid consisting of a two component mixture of 
soft spheres in three dimensions. As the temperature is lowered (or as the density is increased), 
the generalized compressibility drops sharply at the glass transition, with the drop becoming more 
and more abrupt as the measurement time increases. At our longest measurement times, the drop 
occurs approximately at the mode coupling temperature Tc- The drop in the linear generalized 
compressibility occurs at the same temperature as the peak in the specific heat. By examining the 
inherent structure energy as a function of temperature, we find that our results are consistent with 
the kinetic view of the glass transition in which the system falls out of equilibrium. We find no size 
dependence and no evidence for a second order phase transition though this does not exclude the 
possibility of a phase transition below the observed glass transition temperature. We discuss the 
relation between the linear generalized compressibility and the ordinary isothermal compressibility 
as well as the static structure factor. 
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I. INTRODUCTION 

The glass transition is still not well understood de- 
spite extensive study. There have been two main theo- 
retical approaches to the problem: dynamic and thermo- 
dynamic. Theories in the first category view the glass 
transition as a kinetic phenomenon characterized by a 
growing relaxation time and viscosity ^, |j, EL . When 
the relaxation time exceeds the measurement time, par- 
ticle motion appears to be arrested resulting in the glass 
transition. One of the most prominent theories espousing 
this view is the mode coupling theory (MCT) in which 
ideally the relaxation time diverges at a temperature Tc 
above the experimental glass transition temperature ||. 
Interesting and fruitful concepts such as dynamic inho- 
mogeneities [@, |], @] and the influence of the energy land- 
scape on relaxation processes [|| |9| have resulted from 
this approach. The thermodynamic viewpoint attributes 
the glass transition to an underlying phase transition hid- 
den from direct experimental observation by extremely 
long relaxation times Jj|, p|, [To[ [lj] , p"2| ] . In most scenarios 
there is an underlying second order phase transition asso- 
ciated with a growing correlation length which produces 
diverging relaxation times as well as diverging static sus- 
ceptiblities H H H HI, 0, |§- More recently Mezard 
and Parisi p2l |1S|| have argued that the underlying tran- 
sition is actually a random first order transition signaled 
by a jump discontinuity in the specific heat. 

Experimentally the glass transition is characterized by 
both kinetic and thermodynamic features. For example 
in the supercooled liquid kinetic quantities such as the 
viscosity and relaxation time grow rapidly as the temper- 
ature is lowered. When the system falls out of equilib- 
rium below a certain temperature, thermodynamic quan- 
tities exhibit features reflecting the glass transition. For 
example as the system is cooled the specific heat has a 
step-like form and the dielectric constant exhibits a peak 
at a frequency dependent temperature. 



In an effort to better characterize the glass transition 
we introduce a novel probe which we call a generalized 
compressibility (^0|. Unlike the specific heat which mon- 
itors energy fluctuations, this linear compressibility is a 
function of the microscopic structure of the system: it 
depends solely on the positions of the particles and not 
on their previous history. Since we do not need to com- 
pare the system's state at different times, it is not a dy- 
namic or kinetic quantity. Rather it is a thermodynamic 
quantity in the sense that it is purely a function of the 
microstate of the system dictated by its location in phase 
space. The generalized compressibility is easy to compute 
numerically, and it is simpler than the dielectric constant 
which involves both the translation and orientation of 
electric dipoles. In addition it does not suffer from finite 
size effects that can often plague measurements of the 
ordinary compressibility deduced from simulations. The 
generalized compressibility can be calculated in either the 
canonical or grand canonical ensembles. In particular it 
is well defined for a system with fixed volume V and 
particle number N in contrast to the ordinary compress- 
ibility which is defined for a system that has fluctuations 
in N or V . The generalized compressibility should be 
directly measurable experimentally in colloidal suspen- 
sions of polystyrene spheres E| and possibly in other 
systems as well. In this paper we present measurements 
of this quantity in a molecular dynamics simulation of a 
two component system of soft spheres. We find that the 
linear generalized compressibility drops sharply as the 
temperature decreases below the glass transition temper- 
ature T g . The drop becomes more and more abrupt as the 
measurement time increases. This is consistent with the 
structural arrest associated with a kinetic transition in 
which the system falls out of equilibrium. Similar results 
are seen as the density is increased at fixed temperature. 

The paper is organized as follows. Section II de- 
scribes the molecular dynamics simulations. Section III 
describes how the relaxation times and mode coupling 
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Tq are determined. These are useful for setting the time 
and temperature scales. Section IV describes our specific 
heat measurements which show a peak at the glass tran- 
sition. Section V derives the expressions for the linear 
and nonlinear generalized compressibilities, and shows 
our results for these quantities. The linear generalized 
compressibility shows an abrupt drop at the same tem- 
perature and density as the peak in the specific heat. Sec- 
tion VI compares the ordinary isothermal compressibility 
with our linear generalized compressibility and shows the 
advantages of the latter. Section VII gives our results for 
the diffusion constant. Section VIII explains the relation 
between the linear generalized compressibility and the 
static structure factor. Finally we summarize our results 
in Section IX. A brief description of some of these results 
as well as results for a single component fluid that forms 
a crystal was reported earlier 63]. 



II. MOLECULAR DYNAMICS SIMULATION 

We have performed a molecular dynamics simulation 
on a glass forming liquid |^3|, |24j] consisting of a 50:50 
binary mixture of soft spheres in three dimensions. The 
two types of spheres, labelled A and B, differ only in their 
sizes. The interaction between two particles a distance r 
apart is given by V a p{r) = e[(t7 Q/3 /r) 12 + X a p(r)] where 
the interaction length a a [j — (cr a + <Jp)/2 with <tb/&a = 
1.4 (a, j3 = A, B). For numerical efficiency, we set the 
cutoff function X a p(r) = rja afi - A with A = 13/12 12/13 . 
The interaction is cutoff at the minimum of the potential 
Vafiir). Energy and length are measured in units of e 
and a a, respectively. Temperature is given in units of 
e/kg where kg is Boltzmann's constant, and time is in 
units of aA\fm/e where to, the mass of the particles, is 
set to unity. The equations of motion were integrated 
using the leapfrog method Q with a time step of 0.005. 
During each run the average density p a = N/L 3 was 
fixed, and the temperature was kept constant using a 
constraint algorithm p5j| . The volume V = L 3 . N — 
Na + Nb is the total number of particles. The system 
occupies a cube with dimensions (± L/2, ± L/2, ± L/2) 
and periodic boundary conditions. 

We have done sweeps of both temperature and density. 
We fix the parameters so that crystallization is avoided 
upon cooling or when the density is increased. For the 
temperature sweeps, we fix the density at p Q — 0.6. For 
<JbI g a = 1-4, this corresponds to a packing fraction of 
1.04. Having a packing fraction larger than 1 means that 
each particle was interacting with other particles most, if 
not all, of the time. Our measuring procedure is the fol- 
lowing. For runs where we cool the system, we start each 
run at a high temperature (T=1.5) and lower the tem- 
perature in steps of AT = 0.05. At each temperature we 
equilibrate for 10 4 molecular dynamics steps (md steps) 
and then measure the quantities of interest for N T addi- 
tional md steps where N T = 10 5 , 2 x 10 5 , 10 6 , 3 x 10 6 , or 
10 7 . All the particles move at each md step. The results 



are then averaged over up to 40 different initial conditions 
(different initial positions and velocities of the spheres). 
We have done some runs in which we heat the system of 
particles by starting at our lowest temperature T = 0.1 
with a configuration obtained by cooling the system. We 
then increased the temperature in steps of AT = 0.05. 
As before we equilibrate at each temperature for 10 time 
steps and then measure quantities for an additional 10 6 
time steps. 

We have also done some density sweeps in which we 
fix the temperature (T=1.0) and systematically change 
the density. The glass transition occurs as we increase 
the density. Colloidal experiments often study the glass 
transition as a function of density. We start each run at 
a low density (p = 0.4) and increase the density in steps 
of Ap = 0.025. At each density we equilibrate for 10 4 
md steps and then measure the quantities of interest for 
N T additional md steps. 

The glass transition occurs either as the temperature 
is lowered or as the density is raised. It is worth noting 
that temperature and density can be combined into a 
dimensionless parameter T p6fl : 

T = P a eff/ T * °iff =^2nc t npal f} (1) 

a/3 

where u e ff represents an effective diameter for particles 
in the mixture. The concentration of each type of particle 
is given by ua — Na/N and ns = 1 — tia- For our simu- 
lations ha = tib = 0.5. r is the relevant parameter when 
the particles spend most of their time sampling a nonzero 
interparticle potential, i.e., for p -1 / 3 < <J e ff- Thus T is 
particularly useful for interparticle interactions which fall 
off with distance as a power law and do not have a cutoff 
beyond which the interaction is zero. When cooling from 
the liquid phase, the glass transition is known to occur 
around T — 1.45 

We have looked for phase separation of the two types of 
spheres by examining the distribution of large and small 
spheres in the neighborhood of large spheres and in the 
neighborhood of small spheres. We see no evidence for 
phase separation at either high (T = 1.5) or low (T = 
0.15) temperatures at a density of p = 0.6. 

III. RELAXATION TIMES AND MODE 
COUPLING T c 

As points of reference for the time and temperature 
scales, it is useful to find the mode coupling Tc and the 
a relaxation times. We can find the relaxation times us- 
ing the intermediate scattering function F(k, t) which is a 
useful probe of the structural relaxation. It is the spatial 
Fourier transform of the van Hove correlation function 
G{r, t) and the inverse time transform of the dynamic 
structure factor S(k 1 u>). There are two different types 
of intermediate scattering function: the self (incoherent) 
intermediate scattering function F s (k, t) and the full (co- 
herent) intermediate scattering function F(k,t). 
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In a computer simulation, the self (incoherent) part of 
the partial intermediate scattering function F s ^ a {k,t) can 
be calculated directly using |^7| 

where the subscript a refers to the particle type, A or 
B. fi(t) is the position of particle i at time t, and (...) 
refers to an average over different configurations. The 
wave vector k = 2nq/ L where q is a vector of integers. 
For an isotropic system F s ^ a (k,t) depends only on the 

magnitude k = \k\. We will choose k = k max where k max 
is the position of the first maximum of the partial static 
structure factor S a (k). In Figure [l] we show the self in- 
termediate scattering function F Sj B(k,t) versus time at 
temperatures below the caging temperature (T w 0.4). 
The caging temperature is the highest temperature at 
which a plateau is present in the intermediate scattering 
function versus time. The plateau represents the tempo- 
rary localization of a particle by a cage of other particles 
surrounding it. 

Mode coupling theory is applicable in the temperature 
range below the caging temperature and somewhat above 
the mode coupling Tc- We define the relaxation time r s 
by F s (k,T s ) — 1/e. We determine the relaxation times 
for the seven highest temperatures shown in Figure [l| 
and then fit the temperature dependence of r s (T) to the 
mode coupling form t s (T) ~ (T — Tc) -7 to find Tc- For 
the self part of the intermediate scattering function, the 
actual value of t s increases as the magnitude of the wave 
vector decreases Q. However, the value of Tc is inde- 
pendent of k. t s (T) versus temperature and the mode 
coupling fit are shown in Figure |2|. We find the best fit 
with the mode coupling temperature Tc — 0.303 which 
corresponds to T = 1.46. Note that Tc is determined 
from measurements made at temperatures where the sys- 
tem is equilibrated. Also shown in Figure ^ is the fit to 
the Vogel-Fulcher form r s (T) = A exp[B / (T -T v F )] with 
Tvf — 0.21 which corresponds to T = 1.60. In doing the 
Vogel-Fulcher fit, we were able to use a much broader 
range of temperatures since Tvf is much lower than the 
mode coupling Tc- 

The full intermediate scattering function F(k,t) is 
given by {27) 

F a &t) = j^(p~ k Jt)P-kj0)) (3) 

where the Fourier transform of the density pt{£) = 

53^-1 e~ tk ' ri( - t \ The subscript a refers to the particle 
type, A or B. The longest a relaxation time can be de- 
termined from the full intermediate scattering function 
evaluated at k = k max p9| |. We set k max — 27r • 8.3666/L 
(L = 8) which is the location of the first peak in the 
structure factor for type B particles. We define the a 
relaxation time r as the time where F a (k max ,t) decays 
to 1/e. At a temperature T = 0.2895858 which is just 
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time [md steps] 

FIG. 1: The self intermediate scattering function versus time 
for a system with 512 particles and p ~ 0.6. The time is 
given in units of molecular dynamics time steps. From left to 
right, the curves are for temperatures T = 0.381679, 0.373134, 
0.364964, 0.357143, 0.34965, 0.342466, 0.33557, 0.328947, 
0.321543, 0.3021148, 0.289855 respectively. 256 type B parti- 
cles were used and the wave vector k = 2-k ■ 8.3666/L which is 
the location of the first peak in the structure factor for type 
B particles. L = 8. For each curve the system was initialized 
from a configuration at that temperature obtained from par- 
allel tempering which is described in the Appendix. Then the 
simulation was run only at that temperature. The tempera- 
tures were chosen so that the parallel tempering acceptance 
rates were high. The curves at the 7 highest temperatures 
were equilibrated for 1 million md time steps before record- 
ing the configurations used to calculate F s (k, t). Each curve 
of the 7 highest temperature curves is averaged over 24 runs 
except for T = 0.373134 which is averaged over 54 runs. The 
curves for T = 0.328947 and 0.321543 were equilibrated for 2 
million md time steps before recording the configurations used 
to calculate F s (k,t). These two curves were averaged over 11 
runs. The curve for T = 0.3021148 was averaged over 22 runs 
and was equilibrated for 10,000 md time steps before record- 
ing the configurations used to calculate F s (k,t). The curve 
for T = 0.289855 was averaged over 36 runs and equilibrated 
for 50 million md time steps before recording configurations 
used to calculate F 3 (k,t). 



below the mode coupling Tc, we find that Fs(k,t) has 
fallen to 1/e at r = (1.0 ±0.1) x 10 6 md time steps for a 
system with 512 particles of which half are type B. This 
gives us a time scale by which to compare other times 
such as our run times. This value of r shows no signs of 
aging |3(], |3lJ and stays about the same even after 10 s 
time stepsr"At higher temperatures this relaxation time 
is much shorter. 

The runs used to determine the intermediate scatter- 
ing function were done in a slightly different way from 
the other measurements. These runs were performed at 
a given temperature and density for N T md time steps 
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FIG. 2: Relaxation times t s versus temperature. The solid 
line is the mode coupling fit to the form t s = A(T — Tc) 1 
with Tc = 0.303, 7 = 1.735 and A = 47.6. The dashed line is 
the fit to the Vogel-Fulcher form t 3 (T) = Aexp[B/(T — Tvf)] 
with T VF = 0.21, A = 33.3, and B = 0.803. 



with no change in temperature or density. The runs were 
started from a configuration that had been equilibrated 
at that temperature and density using parallel temper- 
ing. The parallel tempering technique is described in the 
appendix. 



IV. SPECIFIC HEAT 

The specific heat is a thermodynamic quantity which 
undergoes a change signaling the glass transition. In ex- 
perimental systems under constant pressure the specific 
heat exhibits a smooth step down as the temperature is 
lowered through the glass transition. In our simulations 
which are done at constant volume, the specific heat has 
a peak at the glass transition. It is a useful check of our 
calculation to see if the peak occurs at the same temper- 
ature (or density) as the drop in the linear generalized 
compressibility. There are two ways to compute the spe- 
cific heat Cv per particle at constant volume V. The first 
is by taking a derivative of the average energy (E) per 
particle with respect to temperature: Cy = d(E)/dT. 
Since we study the system at discrete temperatures, we 
approximate the derivative by a finite difference: 



Cy(T n 



(E(T n )) - (E(T n -i)) 



(4) 



where T n > T n —x for all integers n. The second way to 
calculate the specific heat is from the fluctuations: 



C v = Nk B (3 2 ((E 2 P ) - (E P ) 2 ) 



(5) 



FIG. 3: Specific heat at constant volume as a function of tem- 
perature for binary mixture of 512 particles with a measuring 
time of 3 x 10 6 md steps averaged over 6 runs. p = 0.6 and 
cfb/<7a = 1.4. The specific heat is calculated from energy 
fluctuations and by taking the derivative of the energy with 
respect to temperature. 



our three dimensional simulations the kinetic energy per 
particle is given by 3k B T/2, so it is the fluctuations in 
the potential energy Ep per particle which determine the 
temperature dependence of the specific heat. Thus 



C v = -k B + Nk B (3 2 ((E 2 p) - (Ep) 2 ) 



(6) 



where kg is Boltzmann's constant, f3 is the inverse tem- 
perature, and Ep is the potential energy per particle. In 



In equilibrium these two ways of calculating the specific 
heat should agree. So we compare the results of calculat- 
ing Cy both ways as a check on our calculation and to 
make sure the system has equilibrated in all the basins 
that were visited in the energy landscape. 



A. Specific Heat Versus Temperature 

The specific heat at constant volume exhibits a peak 
at the glass transition as shown in Figure ||. The data 
in this figure is for 512 particles and was averaged over 
6 runs with a measurement time of 3 x 10 6 md steps. 
Notice that there is good agreement between calculat- 
ing the specific heat by taking a derivative of the energy 
with respect to temperature (see eq. @)) and by using 
fluctuations (see eq. @). This implies that the system 
has equilibrated within the basins that it visits in the 
energy landscape. We find similar agreement for other 
run times. At low temperatures the specific heat goes 
to 3k B as expected for classical oscillators while at high 
temperatures Cy approaches 3k B /2 which corresponds 
to an ideal gas. The peak in the specific heat occurs at 
T w 0.3 which corresponds to Y w 1.46. The temperature 
of the peak coincides with the mode coupling Tc = 0.303 
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FIG. 4: Specific heat at constant volume as a function of tem- 
perature for binary mixture of 512 particles with measuring 
times of 10 5 , 2 x 10 5 , 10 6 , 3 x 10 6 , and 10 7 md time steps. The 
number of runs averaged over is indicated in the legend. The 
specific heat is calculated from energy fluctuations. p — 0.6 
and ob I = 1.4. 



FIG. 5: Specific heat at constant volume during heating and 
cooling a binary mixture of 512 particles with a measuring 
time of 10 6 md time steps averaged over 10 runs. p = 0.6 
and Ob jo a = 1.4. 



that we deduced from the intermediate scattering func- 
tion data. Longer run times lead to a sharper peak in the 
specific heat as can be seen in Figure |] which shows the 
specific heat for 512 particles for several different mea- 
suring times. The peaks would presumably be sharper if 
we had used a finer temperature scale. At high temper- 
atures the agreement between the different times is very 
good. Perera and Harrowell Q have found a specific 
heat peak in a two dimensional binary mixture of soft 
spheres. They argue that their peak is an equilibrium 
feature. However, in our case, at temperatures below the 
peak, the system has fallen out of equilibrium and has 
become trapped in a basin in the energy landscape. We 
shall see this later by examining the energy of the inher- 
ent structures (potential energy minima) as a function of 
temperature. Thus the fact that the peak in the specific 
heat occurs at or very close to the mode coupling Tc is a 
result of the relaxation times (see Fig. |^) becoming com- 
parable to and exceeding the simulation run times as the 
temperature drops below Tc- When this happens, the 
system falls out of equilibrium and undergoes a kinetic 
glass transition. 

The specific heat Cp of experimental systems at con- 
stant pressure exhibits a downward step at the glass 
transition during cooling and a peak at slightly higher 
temperatures upon heating [ |33| . As can be seen in Fig- 
ure |, in our warming up simulations, which are done 
at constant volume, the specific heat peak sharpens and 
moves toward higher temperatures compared to the cool- 
ing runs. This is consistent with what is seen in ex- 
periments. The hysteresis is consistent with the system 
falling out of equilibrium and getting stuck in a basin of 



the energy landscape. 

As we mentioned in the introduction, some have sug- 
gested that the glass transition has an underlying second 
order phase transition @ [TJ, [l| [l6[ |TJ, [0|. Unlike 
typical second order phase transitions, there is no exper- 
imental evidence that the specific heat diverges at the 
glass transition. This is consistent with our simulations. 
In simulations one looks for a divergence by examining 
whether the quantity increases systematically with sys- 
tem size. In Figure ^ we plot Cy for systems with 64, 
216, 512 and 1000 particles. As one can see, the specific 
heat does not exhibit any size dependence. However, 
we cannot rule out the possibility that a thermodynamic 
phase transition occurs at temperatures below where we 
fall out of equilibrium. Indeed theories which postulate 
a thermodynamic transition put the transition tempera- 
ture well below the mode coupling Tc- 



B. Specific Heat Versus Density 

In Figure we show the specific heat as a function of 
density. As the density increases, the specific heat rises 
to a peak at pP° ak = 0.8. This corresponds to T = 1.44 
which is in good agreement with the T value of 1.46 that 
we found for the specific heat peak when we varied the 
temperature. Going to higher densities corresponds to 
going to lower temperatures. At densities higher than 
0.8, the system falls out of equilibrium. 
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FIG. 6: Specific heat during cooling for binary mixtures of 
64, 216, 512, and 1000 particles. The measuring time was 
3 x 10 6 md time steps. The specific heat was calculated from 
fluctuations and averaged over the number of runs indicated 
in the legend. Note the lack of size dependence. p — 0.6 and 
cb/cta = 1.4. 



shots of the configurations of the particles at different in- 
stances, scramble the order of the snapshots, and still be 
able to calculate the generalized compressibilities. Aver- 
aging over these snapshots corresponds to ensemble aver- 
aging. In this sense the generalized compressibilities are 
thermodynamic quantities which can be calculated solely 
from the microstates of the system and do not depend on 
the system's dynamics or kinetics. 

We now derive expressions for the linear and nonlinear 
generalized compressibility. To probe the density fluc- 
tuations, we follow the approach of linear response the- 
ory and consider applying an external potential ^f-4>{r) 

which couples to the local density p(r) = ^i=i S(r — fj) 
where denotes the position of the i th particle. p Q is the 
average density. AP has units of pressure and sets the 
magnitude of the perturbation. (f>(r) is a dimensionless 
function of position that must be compatible with the pe- 
riodic boundary conditions imposed on the system, i.e., 
it must be continuous across the boundaries, but is oth- 
erwise arbitrary. This adds to the Hamiltonian H of the 
system a term 



U 



AP 



Po Jv 



d 3 r(f>(r)p(r) 



AP 

Po 



£*(*) 



AP 

Po 



P<P (7) 
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FIG. 7: Specific heat versus density for a binary mixture of 
512 particles with T = 1. The measuring times were 2 x 
10 5 , 10 6 and 3 x 10 6 md time steps. The specific heat was 
calculated from fluctuations and averaged over the number of 
runs indicated in the legend. p = 0.6 and o\b/o"a = 1.4. 



V. GENERALIZED COMPRESSIBILITIES 

As we mentioned in the introduction, the generalized 
compressibilities are thermodynamic probes that are a 
function of the microscopic structure of the system. They 
are solely a function of the positions of the particles and 
do not depend on their histories. So one could take snap- 



where we have defined p^, — J v d 3 r(p(r)p(r) = 4>{fi). 
p$ is the inner product of <j> and p(r), and we can regard it 
as a projection of the density onto a basis function (j>(r), 
i.e., p^ =< p\4> >. ft weights the density fluctuations 
according to their spatial position. The application of 
the external potential will induce an average change 



<W = {P<t>)u 
where the thermal average 

(P<t>)u = -^Tr 



(P4>)u=o 

u is given by 

-0{H+U), 



(8) 



(9) 



The partition function Z = Tre _/3 ( H+C/ ) and j3 is the 
inverse temperature. For small values of AP, this change 
can be calculated using perturbation theory p4| . Up to 
third order in AP, we find 



tip* 



Po 

/3 3 AP 3 



0/C 



P 2 AP 2 

~^pT 



(Pl)c 



6p 3 a 



{pV)c 



where the cumulant averages are 

(pI)c = (pD ~ W) 2 

{P%)c = (4> - 3<p^><4) 
<4>c = (4>-4^><4>-3<4> 2 + 
12( P0 } 2 <4)_ 6 (^)4 

with the thermal average (p^) = (p^)u=o- The third or- 
der cumulant, eq. (p^) , is zero in the liquid phase because 



(10) 

(11) 
(12) 

(13) 
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for every configuration there exists an equivalent config- 
uration with the opposite sign of (p$ — (p$)) and so we 
will not consider this term any further. We can recast 
cq. ( fTo|) as a power series in the perturbation AP: 



$P<f> 
N 

where 



1 



6p a k B T 



1 



Q{pok B Tf 



Xm(AP) 3 (14) 



Xi 



Xnl 



(15) 



In this paper we will focus our attention on the linear 
(Xi) an d nonlinear (xnl) dimensionless generalized com- 
pressibilities defined by the above expressions. 

We now discuss the choice of the function <j>. We con- 
sider applying the potential along the direction p of one 
of the coordinate axes so that <j>{r) — <f>(r^). A natural 
candidate for </>(r M ) is cos(fc M r A ') (no implied sum over 
repeated indices) with k^ = 2Tm/L, where n = 1,2, ... In 
this case, p^ is the k th mode of the cosine transform of 
the density. We will also consider the simpler function 



\/L. The absolute value corresponds to the 



case where all the particles feel a force along the pth di- 
rection pointing towards the origin. It gives results very 
similar to (j){r^) — cos(fc AI r A ') for small k at a fraction of 
the computational cost. (No sum over repeated indices.) 
So our results in this paper correspond to two cases: 



Ewi/ £ 



which is rather like a center of mass, and 



P<t> 



E cos(fc A1 rf) 



(16) 



(17) 



Since the system is isotropic, we compute the compress- 
ibilities for each direction and then average over the di- 
rection p. 

In most of our calculations we work in the canonical 
ensemble where we fix the volume V, the number N of 
particles and the density p . However, it is straightfor- 
ward to generalize our results to the grand canonical en- 
semble where the number of particles is not fixed. We 
simply replace the thermal average defined in equation 
(PI) by 



1 



W)tf = ^$>^Tr 

N 



(H N +U N ) / 



(18) 



where p is the chemical potential, ifyv is the Hamiltonian 
with TV particles, Un is given by eq. (^) for a system 
with N particles, and Z is the grand canonical partition 
function given by 



2 = Yl eMJVTr 



-f3(H N + U N ) 



(19) 



The generalized compressibilities can be defined using 
equations ([ll]) through (JlJ) with the thermal averages 
(Ps) ~ (P4,)u=o defined in the grand canonical ensemble. 
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FIG. 8: Linear generalized compressibility as a function of 
temperature for different measuring times t M : 10 5 (A, 40 
runs), 2 x 10 5 (o, 32 runs), 10 6 (□, 10 runs), 3 x 10 6 (O, 
b runs) and 10 7 (*, 6 runs) md steps. System size is 512 
particles. p = 0.6 and ob/oa = 1.4. xi i s calculated using 
the absolute value of the particles' positions. Inset: T > T 
subset of the same data scaled as described in text. 



A. Results for Linear Generalized Compressibility 

We now turn to our results for the binary glass forming 
liquid. 



1. Xi from Absolute Value of Positions versus Temperature 

We will first discuss the linear generalized compress- 
ibility calculated from the absolute values of the particle 
positions using eqs. ( [l5|) and dlrf). Figure || shows the 
linear generalized compressibility as a function of tem- 
perature for different run times. The compressibility at 
high temperatures is independent of T. In the vicinity 
of the glass transition xi drops. Notice that as the mea- 
suring time tM increases (and hence as the cooling rate 
decreases) , the temperature of the drop decreases and be- 
comes more abrupt. The measuring time can be thought 
of as the number of snapshots at a single temperature 
that we use to calculate the compressibility. The linear 
compressibility is proportional to the width of the distri- 
bution of p$, so the drop in xi corresponds to the sudden 
narrowing of the distribution P(p$). If we regard p^ as 
a generalized center of mass, then the drop in xi signals 
the sudden arrest in the fluctuations of the generalized 
center of mass. In other words, at the glass transition the 
motion of the particles is largely frozen and hence, the 
generalized center of mass does not move around much. 
This is consistent with recent observations of the colloidal 
glass transition in which the size of the clusters of "fast" 
particles drops dramatically at the glass transition |2l|. 



8 



Notice that at longer measuring times, the tempera- 
ture Xd rop at which the generalized linear compressibil- 
ity drops is roughly at the mode coupling temperature 
Tc = 0.303. Let us define Td rop as the temperature at 
which xi has dropped halfway down. For 10 6 md steps, 
Tdrop ~ 0.33; for 3 x 10 6 md steps, T dr0 p « 0.30; and for 
10 7 md steps, Td rop ~ 0.27. Thus we are able to stay in 
equilibrium down to the mode coupling temperature for 
our longer runs. This is what we would expect when we 
compare these run times, which are longer than 1 million 
time steps, to the a relaxation time r which is about 1 
million time steps at T — 0.29 which is just below Tc- 
Thus the fact that the drop in the linear generalized com- 
pressibility occurs at or very close to the mode coupling 
Tc is a result of the relaxation times (see Fig. ||) be- 
coming comparable to and exceeding the simulation run 
times as the temperature drops below Tq. When this 
happens, the system falls out of equilibrium and under- 
goes a kinetic glass transition. 

The behavior exhibited by \i can be quantified us- 
ing a scaling ansatz: xi{^m,T) = g(jx = t M /r(T)), 
where the characteristic time has the Volgel-Fulcher 
form t(T) = exp(A/(T - T Q )). The inset of Figure | 
shows that the data collapse onto a single curve with 
A = 0.75, T = 0.15. (The data could not be fitted using 
t(T) = A{T — To) 1 as suggested by simple mode cou- 
pling theories [||.) Notice that T D lies below the Vogel- 
Fulcher temperature Tvf — 0.21 and the MCT critical 
temperature Tc = 0.303 deduced by fitting the temper- 
ature dependence of the relaxation times. This scaling 
suggests that \i becomes a step function for infinite tM 
and that the drop in the compressibility would become a 
discontinuity at infinitely long times. This abrupt drop 
is consistent with a sudden arrest of the motion of the 
particles in the liquid which is the kinetic view of the 
glass transition. The abrupt drop also appears to be in 
agreement with Mezard and Parisi's proposal that the 
glass transition is a first order phase transition with a 
jump in the specific heat |l2|, ^9[ . Indeed the tempera- 
ture at which xi drops agrees with the temperature of 
the peak in the specific heat shown in Figure § The 
specific heat provides an independent check of the glass 
transition temperature. However, the drop in xi & n d the 
specific heat peak are due to the system falling out of 
equilibrium, and therefore we cannot really tell if there 
is an underlying true thermodynamic transition. 

One way to see that the system is falling out of equilib- 
rium is to plot the inherent structure energy per particle 
H |3^, [56). An inherent structure is a particular system 
configuration whose energy corresponds to the minimum 
of a basin in the energy landscape. The energy landscape 
is a 3N dimensional surface defined by the potential en- 
ergy of the system which is a function of the particles' 
coordinates. During each run we sampled the configura- 
tions found at each temperature. Each configuration lies 
somewhere in a basin and we used the method of conju- 
gate gradients |f37f to find the inherent structure energy 
of that basin. The result is shown in Figure ^ where 
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FIG. 9: Inherent structure energy per particle as a function 
of temperature for a system of 512 particles at different mea- 
suring times. Other parameters are the same as in Figure 



we plot the average inherent structure energies versus 
the temperature of the configuration that was originally 
saved. At high temperatures the inherent structure en- 
ergy e/s per particle is flat as a function of temperature. 
As the system is cooled, eis decreases rather steeply [||. 
The inherent structure energy flattens off at low temper- 
atures where the system has fallen out of equilibrium and 
has become stuck in one basin. For each measuring time 
the temperature below which the generalized linear com- 
pressibility drops corresponds to the temperature below 
which the inherent structure energy flattens off at low 
temperatures. Thus the temperature of the drop in xi 
and the peak in the specific heat corresponds to the tem- 
perature below which the system falls out of equilibrium 
and ceases to explore deeper basins of the energy land- 
scape. 

The behavior of the linear generalized compressibility 
seen in Figure || is similar to that seen in measurements of 
the real part of the frequency dependent dielectric func- 
tion e'{ijj) pSj . In that case as the frequency decreased, 
the temperature of the peak in e'{lS) decreased and the 
drop in s'(lu) below the peak became more abrupt. By 
extrapolating their data to u = 0, Menon and Nagel |ll| 
argued that e'(uj — 0) should diverge at the glass transi- 
tion, signaling a second order phase transition. We have 
looked for evidence of this divergence by examining sam- 
ples of different sizes to see if the linear generalized com- 
pressibility increased systematically with system size. As 
shown in Figure [l^ we find no size dependence and no in- 
dication of a diverging linear generalized compressibility. 
However, as we mentioned earlier, this does not preclude 
the possibility that a thermodynamic phase transition 
occurs below Td r0 p- It simply means that if there is a 
growing correlation length, it is smaller than the size of 
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FIG. 10: Linear compressibility as a lunction of temperature 
for different system sizes: 216, 512 and 1000 particles. The 
measuring time was 3 x 10 6 md steps in all cases. Other 
parameters are the same as in Figure N. 



FIG. 11: Linear generalized compressibility as a function of 
temperature for a binary mixture of 512 particles upon cooling 
and heating. The measuring time was 10 6 md steps in both 
cases. The data was averaged over 10 runs. Other parameters 
are the same as in Figure H. 



our system at T > Td rop . Another possible reason for 
the absence of size dependence may be that if there is 
an underlying thermodynamic phase transition, then its 
order parameter may not couple to the local density p(r) . 

So far we have shown the results of cooling the sys- 
tem. In order to look for hysteretic behavior we have 
done runs in which we heat a system of 512 particles 
by starting at our lowest temperature T = 0.1 with a 
configuration obtained by cooling the system. We then 
increased the temperature in steps of AT = 0.05. As 
before we equilibrate at each temperature for 10 4 time 
steps and then measure quantities for an additional 10 6 
time steps. Our results are shown in Figure 11. Notice 



the slight hysteresis with the rise in \i upon warming 
being at a slightly higher temperature than the drop in 
Xi upon cooling. This hysteresis is consistent with the 
kinetic arrest of motion and with the hysteresis found for 
the specific heat in Figure ||. 



The resulting linear generalized compressibility is qual- 
itatively similar in its temperature dependence to the 
linear compressibility calculated using the absolute val- 
ues of the particles' positions (eq. ([i"6|)). Figure |lj shows 
the linear generalized compressibility versus temperature 
for various values of the wave vector. The data is for a 
binary mixture of 512 particles with a measuring time 
of 10 6 md steps and averaged over 10 runs. Just 
as for the absolute value case, we find that as we in- 
crease the measuring time, the drop in the linear gen- 
eralized compressibility calculated using cosine becomes 
sharper at the glass transition. This is shown in Fig- 
ure [l3] which shows xi as a function of temperature for 
measuring times of 2 x 10 5 , 10 6 , and 3 x 10 6 md steps 
with k = 2n/L, i.e., n = 1. Figure |l4| shows the 
linear generalized susceptibility versus the wave vector k 
in units of 2n/L for various temperatures. Note that the 
dependence is nonmonotonic. 



2. \i from Cosine of Positions versus Temperature 

We now consider calculating the linear generalized 
compressibility from the cosine transform of the density 
using eqs. ( [l5"| ) and (|l7|). So if we apply a cosine potential 
along, say the /i = x direction, then 

p<f> = ^ cosfezj) (20) 

i 

where the wave vector k x — 2nn/L with n = 1, 2, .... The 
wave vectors are compatible with the periodic bound- 
ary conditions of our simulations. Since the system is 
isotropic, we average xi over the x, y, and z directions. 



3- Xi from Absolute Value of Positions versus Density 

In Figure |l5| we plot the generalized linear compress- 
ibility versus density calculated from the absolute value 
of the positions of the particles. In Figure [TJ| we plot 
the generalized linear compressibility versus density cal- 
culated from the cosine of the particles' positions. In 
both cases we see that xi drops with increasing density. 
The drop becomes more abrupt as the measuring time 
increases. This drop is similar to what we see when we 
cool the system at fixed density. The density (p w 0.8, 
r « 1.44) at which the drop occurs agrees with the den- 
sity at which there is a peak in the specific heat as shown 
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FIG. 12: Linear generalized compressibility as a function of 
temperature for a binary mixture of 512 particles for different 
values of the wave vector k = 2irn/L. The measuring time 
was 10 6 md steps in all cases. The data was averaged over 10 
runs. The susceptibility was calculated using eqs. (|TJ|) and 
([[?]). Other parameters are the same as in Figure ^. 



in Fig. |. 



FIG. 13: Linear generalized compressibility as a function of 
temperature for a binary mixture of 512 particles for different 
values of the measuring time. The measuring times are 2 x 10 5 , 
10 6 , and 3 x 10 6 md steps. The data was averaged over the 
number of runs indicated in the legend. The susceptibility 
was calculated using eqs. © and (0). The wave vector 
k = 2-kJL, i.e., n = 1. Other parameters are the same as in 
Figure n2. 



4- Xi from Absolute Value of Positions versus Temperature 
in a Slab Geometry 

So far we have considered systems with a fixed number 
of particles, but as we mentioned earlier in this section, 
we can generalize our results to the grand canonical en- 
semble where the number TV of particles can vary. We 
have examined the generalized linear compressibility xi 
calculated from the absolute value of the particle posi- 
tions using equations ([l5]) and ( |l6| ) for a slab of our sys- 
tem. In other words we have divided a system of 8 3 
particles into 8 slabs of equal thickness perpendicular to 
the x axis. The number of particles in any given slab is 
not fixed. However, in eq. ( p"5| ) we set the average num- 
ber N of particles in each slab equal to the total number 
of particles in the system divided by the number of lay- 
ers, i.e., N = 8 3 /8 = 64. Such a slab geometry mimicks 
experiments on colloidal suspensions of binary mixtures 
in which the focal plane of the camera can essentially see 
only one monolayer of polystyrene balls ]3^| . In figure ^ 
we show the generalized linear compressibility for a slab 
for two different measuring times. Again we see that 
the drop is sharper as the measurement time becomes 
longer. Thus allowing for fluctuations in the number of 
particles does not change the qualitative behavior of \i 
at the glass transition. Comparing Figures || and |l7|, we 
see that the temperature of the drop for the slab and the 
bulk agree. We also notice that the drop is sharper for 
the bulk where presumably the greater number of parti- 



cles results in better statistical averaging. 



5. Xi from Absolute Value of Positions versus Density in a 
Slab Geometry 

Since experiments on colloidal suspensions usually vary 
the density rather than the temperature, we have done 
simulations where we set the temperature T = 1 and vary 
the density. Again we divide our system of TV = 8 3 = 512 
particles into 8 slabs and measure xi m one of those slabs. 

As one can see 



The results are shown in Figure 18 



from the figure, xi drops as the density is increased and 
the drop becomes more abrupt as the measuring time 
lengthens. Comparing Figures [Tj| [l6] and [TJ, we see that 
the drop occurs at the same density (p ks 0.8) as in the 
bulk. 



B. Nonlinear Generalized Compressibility 

We now turn to the case of the nonlinear generalized 
compressibility x-ni given by eq. (fl5|). We are motivated 
by the case of spin glasses where the nonlinear magnetic 
compressibility diverges at the spin glass transition while 
the linear compressibility only has a cusp |l| . There 
have been a few studies of nonlinear response functions in 
structural glasses (IT], [II), but these have not found any 
divergences. Our results are consistent with this con- 
clusion. In particular we find that the nonlinear gener- 
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FIG. 14: Linear generalized compressibility as a function of 
wave vector for a binary mixture of 512 particles for different 
values of the temperature. The temperature is measured in 
md units. The measuring time was 10 6 md steps in all cases. 
The data was averaged over 10 runs. The susceptibility was 
calculated using eqs. ( |l5| ) and (JIt]) . Other parameters are the 
same as in Figure |l2|. 




Density 

FIG. 15: Generalized linear compressibility as a function of 
density for a binary mixture of 512 particles at T = 1. The 
measurement times are 2 x 10 5 , 10 6 and 3 x 10 6 md steps. 
obIoa = 1.4. Data was averaged over the number of runs 
indicated in the legend. \i i s calculated using the absolute 
value of the particles' positions. 



alized compressibility is zero above and below the glass 
transition temperature, though it does show a glitch at 
the glass transition temperature. There is no system- 
atic increase with system size, indicating the absence of 
a divergence at temperatures above the glass transition. 
This does not rule out a divergence below the glass tran- 



FIG. 16: Generalized linear compressibility as a function of 
density for a binary mixture of 512 particles at T — 1. The 
different values of n correspond to different values of the wave 
vector k = 2nn/L. The measurement time is 10 6 md steps in 
all cases. ob/<Ja = 1.4. Data was averaged over 6 runs. The 
susceptibility was calculated from the cosine of the particles' 
positions using eq s. (|l5| ) and (]l7|). Other parameters are the 
same as in Figure |l2|. 
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FIG. 17: Linear generalized compressibility as a function of 
temperature for a monolayer slab in a binary mixture of 512 
particles for different values of the measuring time. The mea- 
suring times are 2 x 10 s and 10 6 md steps. The data was 
averaged over the number of runs indicated in the legend. 
The susceptibility was calculated using eqs. jl5| ) and (fl6"|). 
Other parameters are the same as in Figure B. 
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FIG. 18: Linear generalized compressibility as a function of 
density for a monolayer slab in a binary mixture of 512 parti- 
cles for different values of the measuring time. The measuring 
times are 2 x 10 5 , 10 6 and 3 x 10 6 md steps. The data was av- 
eraged over the number of runs indicated in the legend. The 
susceptibility was calculated using eqs. ( |l5| ) and (|lq). Other 
parameters are the same as in Figure B. 
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FIG. 19: Nonlinear generalized compressibility as a function 
of temperature for binary mixtures. The data for 216 particles 
is from 3.2 x 10 6 md steps obtained by stringing together 16 
runs, each of which involved 2 x 10 md steps. The data for 
512 and 1000 particles is from 6.4 x 10 6 md steps obtained by 



stringing together 32 runs of 2 x 10° 
(Jb/cta = 1-4. 



md steps. p = 0.6 and 



sition temperature where our system has fallen out of 
equilibrium. It also does not rule out a thermodynamic 
transition that does not couple to the local density. Be- 
cause xm is sensitive to the tails of the distribution of 
p^, one must be careful to obtain a good ensemble av- 
erage in the liquid above the glass transition tempera- 
ture. We have done this by doing 16 or 32 runs, each 
involving 200,000 time steps, with different initial condi- 
tions, stringing them together as though they were one 
long run and then taking the appropriate averages. In 
some sense this approach mixes molecular dynamics and 
Monte Carlo; the simulation follows the equations of mo- 
tion for a given amount of time and then "jumps" to 
another configuration which again evolves according to 
the molecular dynamics equations until the next jump. 
We call this approach "global averaging." It produces 
a better ensemble average of < > 2 which enters into 
Xn i in eq. (filf). The resulting Xni is shown in Figure 
|l9| which was calculated from the absolute values of the 
particles' positions using eqs. (JTJ) and (|l6|). 

Xni also took longer to equilibrate than \i ■ By plotting 
Xni versus run time, we found that one had to run at least 
10 6 time steps at T = 1 before \ni appeared to saturate 
(see Figure |2fj|). 



VI. ORDINARY ISOTHERMAL 
COMPRESSIBILITY 

The ordinary isothermal compressibility kt can be re- 
lated to the fluctuations in the number, volume or density 
of the system. In most of our calculations we fix the vol- 
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FIG. 20: Nonlinear generalized compressibility as a function 
of time for binary mixtures of 216 and 512 particles at T = 1. 
The data shown for each system size is for a single run and a 
running average is kept, xm is calculated using eqs. (jl5|) and 
@. pa = 0.6 and a B /a A = 1.4. 



ume, number of particles and density, so that there are no 
such fluctuations and the system has K-p = 0. However 
in the grand canonical ensemble kt is given by 



kt 



1 (N 2 ) - (N) 1 
Pok B T (N) 



(21) 



where p =< N > /V. To relate kt to the linear gener- 
alized compressibility xu we choose a uniform potential 



13 



4>(r) = 1. Then p^ 



J d 3 rp(r) 
1 



N and 



-XI 



(22) 



In principle one can also obtain kt from the k — > 
limit of the static structure factor S'(fe) = (1/N)(p?p_z) 
where p^ — 5Z i=1 exp(— ife ■ r\-) is the Fourier transform 
of the local density p(r). The limit of S(k) for /c — > in 
an isotropic and homogeneous system is |27| 

S(0) = l + P o y ( 5 M-l)d 3 r 

= p k B TK T (23) 

where g(r) is the radial distribution function. Note that 
in a system with fixed volume and particle number, the 
normalization of g(r) leads to S(k — > 0) = 0. This is 
consistent with the fact that such a system has kt = 0. 
Eq. (23|) yields a nonzero value for kt in a system which 
has fluctuations in volume, particle number or density. 
Even in such a compressible system taking the k — > 
limit of S(k) can suffer from finite size effects Eaj be- 
cause the farthest apart that any two particles can be 
along any given coordinate axis is L/2 when there are 
periodic boundary conditions. So at wave vectors k with 
components smaller than 4ir/L, S(k) can have spurious 
results. (For example in our simulations wc found that 
this manifests itself as a slight upturn in S(k) at small k.) 
It is possible to extrapolate to distances larger than L/2 
using various approaches p^ ]. We chose not to use this 
approach to calculate kt since we work in a system with 
fixed N and V. We should note however that simulations 
fi3f working in the NVT ensemble with fixed N, V, and 
T, have successfully used eq. ( |23| ) to find the isothermal 
compressibility. We can resolve this with the fact that 
S(k — > 0) = in the NVT ensemble by noting that for 
values of k > 4:w/L, S(k) should give the same value in 
the NVT ensemble as in the grand canonical ensemble. 
So if L is large enough, fitting S(k) to the small k form 
S(k) = S(k — > 0) + Ak 2 , where A is a constant, should 
yield the correct value of kt as long as k > 4ir/L. 

Rather than using eq. (p3|), we calculated kt by mon- 
itoring a small subvolume inside of our system and keep- 
ing track of the fluctuations in the number of particles 
in the subvolume. Let us define a dimensionless ordinary 
isothermal compressibility Kt by 



1 



where 



kt 



K T = 



pk B T 



(N 2 ) - (N) 2 
W) 



(24) 



(25) 



In order to calculate Kt, we have monitored a subvolume 
that had on average 25% of the total number of particles. 
Essentially we drew an imaginary boundary in the mid- 
dle of our system that enclosed 25% of the total volume 
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FIG. 21: Dimensionless ordinary isothermal compressibility 
Kt as a function of temperature for a subset of a system of 
512 particles for different values of the measuring time. On 
average the subset had 128 particles in it. The measuring 
times are 2 x 10 5 and 3 x 10 6 md steps. The data was av- 
eraged over the number of runs indicated in the legend. The 
dimensionless compressibility was calculated using eq. (pBj). 
Other parameters are the same as in Figure H. 



and counted the number of particles in this subvolume 
as a function of time. By monitoring the fluctuations 
in the number of particles in this subvolume, we could 
calculate Kt- The results for a subvolume which had 
on average 128 particles out of a total of 512 particles 
are shown in Figure ||l] where Kt is plotted versus tem- 
perature. We see that it has the same basic shape as 
the linear generalized compressibililty with a drop at the 
same temperature as Xl- As with xi, the drop becomes 
sharper with increasing measuring time. 

While Kt shows behavior similar to the linear general- 
ized compressibility xi as a function of temperature for a 
given size, it is unlike xi m that it suffers from finite size 
effects. We have demonstrated this by making measure- 
ments on systems with a total of 64, 216, 512, and 1000 
particles. The measurements were made by counting the 
number of particles in a subvolume that was 25% of the 
total volume. Such a small subvolume has a large surface 
to volume ratio which produces large finite size effects. 
To understand this, we note that in such a small subvol- 
ume, a significant number of the particles are very close 
to the boundary of the subvolume. Fluctuations in the 
positions of these particles moves them in and out of the 
subvolume, producing large fluctuations in the number of 
particles in the subvolume. The smaller the system, the 
bigger this effect is. This produces large finite size effects 
in Kt even at high temperatures where the system easily 
equilibrates. This can be seen in Figure |22|. One can see 
that Kt decreases with increasing system size at high 
temperatures above the drop in the compressibility. One 
of the advantages of the linear generalized compressibility 



14 



0.3 



= 0.2 - 



00 

oo 

CD 
Q. 

E 
o 
O 

"cs 

CD 



0.1 - 







V-u-l- 1 Lj-fH i—i i—i — _ i—i i— i |— ljn — ,— . j— i i-tJI ,— , r 1 


/7\ i — i i— • i_i lj-^ ftrn-i_ i l-i — lj u 1=1 Lzn- 1 — trtj 


/ T G— 


— G N=64, 3 x 1 6 md steps (1 runs) " 


/ / 13 — 


— a N=216, 3 x 10 6 md steps (5 runs) 




— N=512, 3 x 10 6 md steps (5 runs) 


/ A— 


—A N=1 000, 3 x 1 6 md steps (7 runs) ' 



10" 



0.5 1 

Temperature 



FIG. 22: Dimensionless ordinary isothermal compressibility 
Kt as a function of temperature for systems with a total of 
N = 64, 216, 512 and 1000 particles. The isothermal com- 
pressibility was obtained by monitoring a subvolume that had 
on average 25% of the particles in it. Notice the size depen- 
dence at high temperatures. The measuring time is 3 x 10 6 
md steps. The data was averaged over the number of runs in- 
dicated in the legend. The dimensionless compressibility was 
calculated using eq. Other parameters are the same as 

in Figure H. 



is the absence of these finite size effects. In fact the lin- 
ear generalized compressibility shows no size dependence 
at tem per atures above the observed glass transition (see 
Figure 



VII. DIFFUSION CONSTANT 

The diffusion of the particles reflects the kinetics of the 
system and becomes very small below the glass transition 
temperature. We calculate the diffusion constant D using 
the Einstein relation 



D 



lim 



N 



(26) 



i=l 



where are true displacements of the «th particle. Since 
we are using periodic boundary conditions, if the parti- 
cle has crossed the box several times, then this must be 
included in r,. 
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FIG. 23: Diffusion constant as a function of temperature for 
a binary mixture of particles plotted on a logarithmic scale. 
Measurement times for 512 particles are 10 5 , 2 x 10 5 , 10 6 , 
3x 10 6 , and 10 7 md steps. Measurement time for 216, 512, and 
1000 particles is 3 x 10 6 md steps. Number of runs averaged 
over is given in the legend. p = 0.6 and (tb/ga = 1.4. 



over the entire temperature range. The curves corre- 
sponding to different cooling rates begin to separate as 
the system falls out of equilibrium at the glass transition 
temperature where the specific heat peaks and where the 
linear generalized compressibility drops abruptly. Figure 
p3| also shows the diffusion constant for binary mixtures 
of several different sizes. Notice that there is no apparent 
size effect. 



B. Diffusion Constant Versus Density 

In Figure |i] we show the diffusion constant as a func- 
tion of density. We see that the diffusion decreases 
smoothly as the density p increases. The curves cor- 
responding to different measurement times begin to sep- 
arate at the glass transition density where the specific 
heat peaks and where the linear generalized compress- 
ibility drops abruptly. 



VIII. RELATION BETWEEN X i AND S(k) 



A. Diffusion Constant Versus Temperature 

As the system is cooled through the glass transition, 
the diffusion constant calculated using equation (26) be- 
comes very small. This is shown in Figure ^3| where the 
diffusion constant for 512 particles is plotted on a log- 
arithmic scale. The diffusion constant varies smoothly 



When the system has translational invariance, we can 
relate the linear generalized compressibility xi to the 
static structure factor S(k) which is measured in experi- 
ments such as neutron scattering. S(k) is also used as an 
input for mode coupling theories ji]] and it is generally 
assumed that S(k) does not show any essential variations 
near the glass transition as the temperature or density is 
varied. As we shall see, our calculation agrees with this. 
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FIG. 24: Logarithmic plot of the diffusion constant as a func- 
tion of density for a binary mixture of 512 particles at T — 1. 
The measurement times are 2 x 10 5 , 10 6 and 3 x 10 6 md 
steps, ob I <JA = 1.4. Data was averaged over the number of 
runs indicated in the legend. 



By definition, the static structure factor S(k) — 
(l/N){pgp_g). To relate S(k) to Xh we nrf =t n °te that 
S(k) is the Fourier transform of the static density-density 
autocorrelation function G(f), i.e., S{k) — J exp(— ik ■ 
r)G(r)d 3 r, where G(r) = (l/N) j (p(f' + r)p(r'))d 3 r' . 
In the supercooled liquid above the glass transition, the 
system has translational invariance, and we can write 



(p(r)p(f')) = —g{r-r') 



(27) 



where V is volume and g is a function of the difference 
(r — r '). In this case G(r) = g(r)/N. In xi we meet 



(pi) = J d 3 rd 3 r^(f)^r'){p(r)p(r')) 
If there is translational invariance, 

1 I" d 3 h 

-Mk)^-k)g{-k) 



(28) 
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V J (2tt) 3 
d 3 k 



d 3 rd 3 r'c/)(r)<j)(r')g(r- r') 
d 3 k 



Po 



(2tt) 3 



\m?s(k) 



(29) 



where 4>(k) is the Fourier transform of 4>{r). Converting 
the integral V f d 3 k/(2i:) 3 to a sum and using eq. 
(|l5|), we obtain 



Xi 



N 



N 



J2(\m\ 2 s(k) 



(30) 



As an example, let us choose 4>(r) — cos(k x x) with 
the proviso that k x ^ 0. (When k = 0, xi = since 
the potential is uniform and there are no fluctuations 
allowed with fixed N and V.) With this choice of <f){r), 
p$ = Y^i C0S (kxXi)- Then one can show by explicitly 
calculating 4>(k) and by using eqs. (^9|) and (||(]) that 



X i(k) = 3S(k) - |p[Re<^}] 2 
where translational invariance allows us to write 



(31) 



S(k) = (2/N)(pl) 

= (2/N)J2[( c °s(kxX i )co8(k x x j ))} (32) 

Note that the value of k that we use to probe the glass 
transition is typically of order k^ = 2it/L which is much 
smaller than the value of fc pea k at which the first peak 
in S(k) appears. fc pC ak ~ ^w/aA 3> fci where a a is the 
diameter of type A spheres in the binary mixture. 

We have numerically calculated [5(fc)] run for our bi- 
nary mixture using equation ( |32| ) with <j>{r) = cos(/c /J r A1 ). 
(No sum over repeated indices, p = x, y, or z.) [...] run is 
an average over all the runs and over k being parallel to x, 
y, and z. The result is shown in Figure ^5| and one can see 
that [S l (fe)] run does not vary much through the glass tran- 
sition which is consistent with what is assumed in mode 
coupling theory. The figure also shows [Re < pk >] rU n 
where < ... > is a thermal average over a single run. If 
[S'(fc)] run and [Re < pk >] rU n do not vary much through 
the glass transition, how can the difference between the 
two terms in ( |3l"| ) decrease and produce a drop in xz? To 
answer this, note that there are two inequivalent ways in 
which one can calculate Xh So far we have calculated Xl 
for each run and then averaged over the different runs. 
This approach is what we used in Figs. and results 

in a sharp drop in the linear generalized compressibility 
at the glass transition. Let us call this a run-by-run 
average for which we can write: 



G 



(33) 



The drop in xi comes about because the width of the 
distribution of Re(p<p) becomes much smaller below the 
transition. At low temperatures structural arrest hinders 
the exploration of phase space and reduces the fluctua- 
tions in He(p^). 

The other way to calculate the generalized linear com- 
pressibility is with global averaging in which we string 
together a series of separate runs, treat it as one giant 
run, and then do the averaging required to calculate the 
generalized linear compressibility xf lobal - 

xf l0bal = ^{[^)] mn -([W)] run ) 2 } (34) 



Careful inspection of eqs. ( |33| ) and (|34|) reveals that the 
difference is in whether we square the thermal average 
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FIG. 25: Linear generalized compressibility Xh static struc- 
ture factor [S(A;)] run , [Re < pk >]run, linear generalized com- 
pressibility xi averaged run by run, and the globally averaged 
linear generalized compressibility global versus temperature 
for a binary mixture of 512 particles. The measurement time 
was 10 6 md steps for each temperature. The data was aver- 
aged over 10 runs and over k — (2n/L, 0, 0), k = (0, 2n/L, 0), 
and k = (0, 0, 2%/L). p a = 0.6 and a B la a = 1-4. 

(ptj,) and then average over runs to obtain [(p<?i) 2 ] run in 
Xh or average (p^) over runs and then square it to obtain 
([W)] m n) m xf lobal - If wc now return to Figure |25| and 
take the difference of [S'(fc)] run and ([Re < pk >] rU n) 2 , we 
obtain the global average: 

x f obal = 3{S(k)] IUU - A ([Re < Pk >] run ) 2 (35) 

The result of both types of averaging is shown in Figures 
p5| and |2(| Notice that xf lobal does not exhibit a drop 
with decreasing temperature while xi does. To under- 
stand why there is no drop in xf'° ba \ note that by com- 
bining several different runs, very different configurations 
are sampled which produces much larger fluctuations in 
the generalized center of mass at low temperatures com- 
pared to xi- As a result xf'° ba \ which is a measure of the 
size of these fluctuations, does not have an abrupt drop. 

IX. SUMMARY 

To summarize, we have introduced a new quantity 
called the generalized compressibility which depends 
solely on the positions of the particles and not on their 
histories. The generalized compressibility can easily be 
calculated in the canonical (e.g., NVT) and grand canon- 
ical ensembles. In particular it is well defined in a system 
which has particle number and volume fixed. In addi- 
tion it does not suffer from the hnite size effects often 



FIG. 26: Linear generalized compressibility as a function of 
temperature for binary mixtures of 216, 512 and 1000 parti- 
cles. The filled symbols correspond to calculating xi f° r each 
run with a measuring time of 2 x 10 5 md steps and then aver- 
aging over the number of runs indicated in the legend, while 
the open symbols correspond to global averages (stringing all 
these runs together to get one big "run" for a given system 
size). So the global average for 216 particles uses 3.2 x 10 6 
md steps while for 512 and 1000 particles, 6.4 x 10 6 md steps 
were used. Xi ls calculated using eqs. © and @. Po = 0.6 
and (J b I o a = 1-4. 



encountered in calculating the ordinary compressibility. 
The linear generalized compressibility drops abruptly at 
the observed glass transition due to the kinetic arrest 
of motion. This makes it an good quantity to calculate 
or measure in order to find the observed glass transition 
as a function of density or temperature. The general- 
ized compressibility can be experimentally measured in 
several ways. It can be directly measured in colloidal ex- 
periments which monitor the positions of the particles. 
Measurements of the width of the distribution of p^, the 
spatial Fourier transform of the density, would also yield 
the linear generalized compressibility. 

We thank Sharon Glotzer, Walter Kob, Andrea Liu 
and Francesco Sciortino for helpful discussions. This 
work was supported in part by DOE grant DE-FG03- 
00ER45843 as well as by CULAR funds provided by the 
University of California for the conduct of discretionary 
research by Los Alamos National Laboratory. 



X. APPENDIX: PARALLEL TEMPERING 

In calculating the intermediate scattering function at a 
given temperature, we initialized the run using a config- 
uration at that temperature generated by parallel tem- 
pering. In this appendix we briefly describe the parallel 
tempering method. 
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We implement parallel tempering (PT) Q ||, |6[ Q 
by choosing the temperatures at which we wish to have 
measurements made. We then run molecular dynamics 
simulations in parallel at these temperatures using a tem- 
perature constraint algorithm p5| to keep the tempera- 
ture of each simulation constant. At 100 time step in- 
tervals we attempt to switch the configurations of two 
neighboring temperatures using a Metropolis test which 
ensures that the energies of the configurations sampled 
at any given temperature have a Boltzmann distribution. 
Let Pi and /?2 be two neigboring inverse temperatures, 
and let U\ and U2 be the corresponding potential energies 
of the configurations at these temperatures at a time step 
just before the possible swap. If A = (/3i — /?2)(C/2 — U%), 
then the switch is accepted with probability unity if 
A < and with probability exp(-A) if A > 0. The 
temperatures are chosen so that the acceptance ratio is 



between 30% and 75%. At the temperatures in the vicin- 
ity of the mode coupling Tq, the acceptance ratio was 
typically above 75% for L=6 and above 60% for L=8. 
After a swap is accepted, the velocities of the particles in 
each configuration are rescaled to suit their new tempera- 
ture. Each configuration is then evolved using molecular 
dynamics for another 100 time steps. Switching config- 
urations allows a given simulation to do a random walk 
in temperature space in which it visits both low temper- 
atures and high temperatures. This helps to prevent a 
simulation from becoming trapped in a valley of the en- 
ergy landscape at low temperatures. Typically we equili- 
brate for 2 million time steps and then do measurements 
for an additional 4 million time steps. 

^Present address: Internap, Seattle, WA 98101. 
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